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ABSTRACT 

We present a more accurate numerical scheme for the calculation of diffusive shock accel- 
eration of cosmic rays using Stochastic Differential Equations. The accuracy of this scheme 
is demonstrated using a simple analytical flow profile that contains a shock of finite width 
and a varying diffusivity of the cosmic rays, where the diffusivity decreases across the shock. 
We compare the results for the slope of the momentum distribution with those obtained from 
a perturbation analysis valid for finite but small shock width. These calculations show that 
this scheme, although computationally more expensive, provides a significantly better perfor- 
mance than the Cauchy-Euler type schemes that were proposed earlier in the case where steep 
gradients in the cosmic ray diffusivity occur For constant diffusivity the proposed scheme 
gives similar results as the Cauchy-Euler scheme. 

Key words: Methods: numerical - Acceleration of particles - Shock waves - Diffusion. 
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1 INTRODUCTION 

One of the leading candidates for the mechanism responsible for 
the acceleration of Galactic cosmic rays is diffusive shock accel- 
eration (DSA). Proposed originally by Krimsky (1977), Axford, 
Leer & Skadron (1977), Bell (1978a,b) and Blandford & Ostriker 
(1978), the theory is now well-understood in the case of non- 
relativistic shocks (shock speed Vj *k c), see for instance the re- 
views of Drury (1983), Blandford & Eichler (1987), Achterberg 
(2001) and Malkov & Drury (2001). In the test-particle limit in 
a steady flow, where the accelerated particles have no influence on 
the flow, DSA at a strong, infinitely thin hydrodynamic shock yields 
a power-law distribution in momentum for the accelerated particles 
with a spectral index q = -dlnM(p)/dlnp ^ 2, where p is the 
particle momentum and n(p) the momentum distribution. Such a 
power law is inferred for cosmic rays at the source between sev- 
eral GeV/nucleon and ~ 100 TeV per nucleon, after a correction 
of the spectrum observed at Earth for the effects of propagation in- 
side (and escape from) the Galaxy. In time-dependent flows with a 
complex flow geometry, in flows that contain multiple shocks and 
in particular when the cosmic rays through their pressure signif- 
icantly decelerate the pre-shock flow, analytical results are more 
difficult to obtain. In those cases one often has to resort to numeri- 
cal methods to solve the basic equations. 

When simulating DSA one essentially tries to determine the 
cosmic ray density N(x , p , t)m reduced phase space (x , p = \p\) 
as a function of time. Two main approaches are possible: one either 
solves the Vlasov equation for the distribution function N{x , p ,t) 
directly, or one uses an algorithm that constructs N(x , p , t) from 
particle positions in phase space that have been obtained by direct 
simulation of representative particle orbits. The first approach re- 
quires the solution of a partial differential equation in (reduced) 
phase space, which is generally computationally expensive as it re- 



quires matrix inversion in schemes such as the Crank-Nicholson 
method (e.g. Potter, 1973) that is needed in order to stably solve 
a diffusion-advection type equation with sufficient accuracy. The 
advantage of course is that one obtains the distribution function 
directly. The second approach is simple to implement in arbitrary 
geometries and can use very accurate integration schemes as one in 
principle solves an ordinary differential equation. Its disadvantage 
is that one must take measures such as particle splitting (see below) 
in order to minimize the effect of Poisson noise that is unavoidable 
as one uses a finite number of particles to construct the distribution 
function. 

This paper is concerned with the second method, which in- 
volves the solution of stochastic differential equations (SDEs). It 
was first proposed in this context by Achterberg & Kriills (1992). 
Recent applications of this method include Marcowith & Casse 
(2010) and Schure et al, (2010). 

We describe a relatively simple stochastic predictor-corrector 
scheme for diffusive shock acceleration. It works well in the pres- 
ence of strong gradients in the diffusivity of the particles. In that 
case this scheme is considerably more accurate than the simple 
Cauchy-Euler scheme, which fails to return accurate results for the 
spectral slope q. The slope returned by the Cauchy-Euler scheme 
is consistently too steep, which indicates that the shock transition 
(where the acceleration takes place) is not sampled accurately. If 
the diffusivity gradient is small or vanishes, both schemes give al- 
most identical results. 

We outline the different numerical approaches in Section 2, 
discuss the need for a better scheme in Section 3, introduce a 
predictor-corrector type scheme in Section 4 and evaluate its per- 
formance in Section 5. Conclusions are found in Section 6. 
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2 NUMERICAL APPROACHES TO DSA USING 
STOCHASTIC DIFFERENTIAL EQUATIONS 

The aim of this paper is to construct the energy distribution of par- 
ticles (cosmic rays) that are accelerated in a prescribed flow con- 
taining a shock, thereby limiting ourselves to the test-particle case. 
The particles undergo spatial diffusion with respect to the flow as a 
result of pitch angle scattering on magnetic fluctuations, and gain 
(or lose) energy due to compression (expansion) of the fluid. Other 
processes, such as radiation losses and stochastic acceleration by 
plasma waves, can be added in a simple manner but will be ne- 
glected here. 

We consider the simplest case of a one-dimensional, steady 
flow along the x-axis with fluid velocity V{x). Let 

AN 

N{x,y,t)^—— (1) 
dx ay 

be the cosmic ray distribution function, where y is a logarithmic 
momentum variable, defined formally as 



<dW,) = 



At . 



(6) 



y = ln(p/rac) 



(2) 



Here m is the particle rest mass and c is the velocity of light. The 
particles are coupled to the fluid by frequent pitch angle scattering 
by scattering centers (magnetic field fluctuations due to hydromag- 
netic waves) that are themselves tied to the flow. In the simplest 
case, where these field fluctuations are advected passively by the 
flow with the local flow velocity V(x), the cosmic ray distribution 
N(x , y , t) satisfies the equation (cf. Skilling, 1975): 



dN d I ^ . .dN\ 1 IAV\ ON 
_,_(v(,).,_0(x)-) = -(-)- 



(3) 



Here N(x , y , t) is the distribution (averaged over gyration phase 
in the laboratory frame, the shock rest frame in the application dis- 
cussed below. D{x) is the position-dependent spatial diffusion coef- 
ficient. For simplicity we assume this quantity to be independent of 
p (and y), but this is not important for what follows. The limitation 
to a steady flow is also not fundamental: the same techniques can 
be used to follow cosmic rays in a time- varying flow. 

It is known that the solutions to this Fokker-Planck type equa- 
tion, which essentially expresses the propagation of particles in the 
two-dimensional phase space (x , y), can be constructed by follow- 
ing particles that satisfy a stochastic differential equation (SDE) 
(e.g. Gardiner, 1983, Ch.5; 0ksendal, 1992, Ch. 3). 

2.1 Spatial transport 

Spatial transport in this case is described by a SDE of the type 
(Achterberg & Kriills, 1992): 



The brackets represent an average over many statistically indepen- 
dent realizations of this Wiener process. 

In practice, this average is achieved by simulating a large num- 
ber of particles with statistically independent orbits. The distribu- 
tion of these particles in phase space approaches the solution of the 
corresponding Fokker-Planck equation, provided one uses a suffi- 
cient number of particles in order to minimize the effects of Poisson 
noise. 

In the context of Diffusive Shock Acceleration the typical 
length scale in the cosmic ray precursor ahead of the shock is the 
diffusion length 



DIV . 



(7) 



This implies that the drift velocity within the precursor is of order 



Vdrift — 



AD(x) AD I AD 



Ax 



D 



(8) 



Taking AD ^ D and V - V^, with Vs the shock speed, one finds that 
Vdrift - ^s- Formally problems arise if V^nft - v with v the particle 
velocity. However, in that case the diffusion approximation that is 
the basis of this method fails. The diffusion approximation presup- 
poses that the anisotropy in the exact momentum distribution of the 
accelerated particles, which is of order K/v, is small. In relativistic 
shocks, where Vs ^ c, this is never the case. When Vdrift ^ c for 
relativistic particles) one has to simulate the pitch angle scattering 
of particles directly, as is routinely done in simulations of particle 
acceleration at relativistic shocks (e.g. Bednarz & Ostrowski, 1998; 
Achterberg et al, 2001). 

Within the shock transition itself one expects a change in dif- 
fusivity AD ^ D over the shock width Ls so that 



D 



'Irafp 

sir 



(9) 



Here /i,^fp is the scattering mean free path so that the diffusion co- 
efficient equals D = V/i^ifp/S. The diffusion approximation requires 
Vdiift V so formally the method fails if Lj < /Imfp. However, 
in that case one may as well approximate the shock as a discon- 
tinuity with a stepwise jump in the diffusivity and flow speed. A 
scaling method for dealing with such sudden step-like jumps in 
the diffusion coefficient and velocity has been devised by Zhang 
(2000). It can be applied in this case so that one solves an equiva- 
lent advection-diffusion problem where no drift term due to the dif- 
fusivity jump at the shock occurs, and the method outlined below 
is valid as long as the diffusion approximation applies to particles 
ahead of and behind the shock front. 



dx = U{x) At + V2D(x) diy, . 



(4) 



In Eqn. (O the quantity 
dD(x) 



t/(x) = V{x) + 



Ax 



(5) 



is an effective advection velocity that includes a drift term due to 
the dynamical friction that results from a spatial gradient in the 
diffusivity. 

The first term in SDE l|4j is a deterministic term oc df that 
describes the average flow of particles while the second stochas- 
tic term represents the spatial diffusion. That last term involves a 
infinitesimal Wiener process dW,, which satisfies 



2.2 Change in particle momentum 

Transport in momentum (in the absence of radiation losses or 
second-order Fermi acceleration by waves) follows from 



An 1 lAV\ 



(10) 



This gives the momentum changes in response to compressions or 
rarefactions in the flow, with (jj{x) the local acceleration rate. The 
effect of radiation losses and of second-order Fermi acceleration 
can be added in a simple manner. 
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2.3 The Cauchy-Euler scheme 



3 THE NEED FOR A BETTER SCHEME 



The simplest numerical scheme that one can use to simulate this 
diffusion-advection process is the explicit Cauchy-Euler scheme 
(CES), see for instance Achterberg & Kriills, 1992. This scheme 
advances x and y in time (time step: A?) with increments in position 
and log momentum Ax and Ay given by: 



Ax = U(x) At + ^j2D(x)At ^, = Axs + Axjiff 



Ay = (jj{x) At . 



(11) 



Here ^, is a normally distributed Wiener process with zero average 
and unit dispersion, in the notation of Kloeden & Platen (1992): 



6 N(0, 1) . 

The quantity 
Axjiff = yj2D(x)At 



(12) 



(13) 



is the rms diffusive step, given the time increment At. The statisti- 
cally sharp (non-stochastic) step is (see Eqn. |5} 



dD 

Ax, = VAt+ — At = Axadv + Axd„f, 
dx 



(14) 



It consists of the advective step Axajv = VAt due to the plasma flow 
and the drift term Axjiif, oc dD/dx due to gradients in the diffusivity. 
This scheme has the virtue of simplicity, and for small Axd,ift it 
can accurately describe cosmic ray acceleration near shocks with a 
thickness provided the time step is chosen in such a way that 



AXadv <C Ls <C AXdi 



(15) 



(Achterberg & Kriills, 1992), who only tested the CES for uniform 
diffusivity so that Axd,.|ft = 0. 

The first condition is necessary to resolve the shock transi- 
tion, and is essentially a demand of sufficient accuracy. The second 
condition ensures that a test particle, while crossing the shock with 
average velocity U and before escaping into the downstream flow, 
crosses the shock many times in diffusive jumps that are typically a 
few times the shock thickness. The simulated particles essentially 
mimic the behavior of a real cosmic ray undergoing acceleration 
near a strong shock, such as a supernova blast wave. 

The whole procedure can be thought of as a replacement of the 
real diffusion (with a microscopic step size) by an equivalent diffu- 
sion process with the same diffusion coefficient but with a macro- 
scopic step size. The last condition then ensures that the acceler- 
ation rate ai(x) due to the compression inside the shock transition 
is sampled frequently and in a stochastic manner. In this way the 
simulation produces the correct spectrum of accelerated particles. 

Marcowith & Kirk (1999) have proposed a partially implicit 
scheme that uses the new particle position x + Ax to evaluate Ay. 
They showed in the simple case of a linear shock transition and 
a constant diffusivity D (so that once again Axjiif, = 0) that this 
approach leads to good results, while the condition (115) can be re- 
laxed to Axadv ^ AX[iiff . Their approach opens the possibility of us- 
ing larger time steps and considering zero-thickness shocks, where 
the velocity change is modeled by a step function. In our approach 
outlined below we will follow their method for the calculation of 
the momentum change in terms of Ay, but will need to keep an 
explicit algorithm for the integration of the SDE for Ax. 



During numerical experiments, where acceleration of test particles 
is calculated in a flow that is obtained using a numerical MHD 
code, it came to our attention that the simple Cauchy-Euler scheme 
does not yield satisfactory results if the drift due to gradients in the 
diffusivity becomes too large, in particular when lAxjiif,! » lAxajvl- 
If the scale length of the variation in the diffusion coefficient is 
Li = \(l/D)(dD/dx)\-^ one has 



Axd, 



Ax^dv 



D 



T7 



(16) 



Here L^ff = D/V is the characteristic length of DSA, the diffusion 
length introduced above. We have found that the CES fails to give 
the correct particle distribution when 



D 



(17) 



In particular (as illustrated below) the momentum distribution of 
the accelerated particles produced by the CES for reasonable val- 
ues of the time step At is too steep. This implies that the average 
momentum gain per shock crossing as calculated using the CES is 
too low and the acceleration rate w(x) is apparently not sampled ac- 
curately by the particles following the simulated orbits. We reiterate 
that Achterberg & Kriills (1992) and Marcowith & Kirk (1999) did 
not consider this case in their calculations: they assumed a constant 
diffusivity D. 

The case of a diffusivity that varies on a scale comparable with 
the shock transition is astrophysically important. For instance: it is 
often assumed that the scattering of cosmic rays near an accelerat- 
ing shock is due to saturated Alfven wave turbulence where the dif- 
fusivity scales as D oc B the case of Bohm diffusion. In a shock in 
an infinitely conducting plasma the MHD shock conditions imply 
that B increases across the shock by a factor 



re 



yjcos- 6b + sin^ 6^ 



(18) 



The parameter r = pa/Pi is the shock compression ratio, and 
9b = cos"'(« • B]) is the inclination angle between the upstream 
magnetic field Bi and the normal h to the shock surface. De- 
pending on the orientation of the upstream magnetic field one has 
I < rB < r. Additional changes in the diffusivity can arise through 
the reflection and transmission of MHD waves at the shock. 

Here (and in what follows) we will use the subscripts 1 (2) to 
describe the value of quantities ahead of (behind) the shock transi- 
tion. A similar case is obtained if the cosmic rays lead to field am- 
plification through the Bell-Lucek instability (Bell & Lucek, 2001): 
there the field is amplified on the diffusive scale Ljiff ~ D/V, in ad- 
dition to the amplification by compression at the shock. 

The use of a fully implicit scheme is not feasible for our ap- 
plication, as the velocity field V(x, t) and the magnetic field fi(x, f) 
are not analytic functions but are determined by a MHD code. We 
therefore need a scheme that is explicit for the advance of the po- 
sition X, that is: it uses the variables at time t and old position x to 
calculate the change Ax and the new position x+ Ax. Such a scheme 
should be sufficiently accurate, numerically stable, and should be 
able to deal with the drift induced by strong gradients in the diffu- 
sivity. 

For the present application it is also important to minimize 
the number of calculations of the spatial derivatives (or avoid them 
altogether) of the velocity field and the diffusivity. Such derivatives 
tend to be noisy when they are determined from the raw output of a 
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MHD code. In the case of the calculation of the momentum change 
this can be achieved by using the method proposed by Marcowith 
& Kirk (1999). In terms of y = ln(p/mc) one replaces the second 
equation of M\\ by: 



Ay: 



[V(x + Aj) - V(x) ] 



;Af . 



(19) 



assuming a steady flow for simplicity. The position change Ax is 
determined in the way outlined below. This algorithm essentially 
uses the spatial average of the acceleration rate along the orbit of a 
simulated particle. For a steady flow: 



Ax j. 



Ax 



I dV(x) 
'3 dx 



V(x + Ax) - V(x) 
3aI 



(20) 



For time- varying flows (where dV/dx 3 V/5x), and in the un- 
likely case that the advective step and the diffusive step cancel each 
other (so that Ax = 0) this prescription can lead to singular be- 
haviour. However, this is easily caught in a numerical scheme and 
correctly dealt with by putting Ay = in that case. 



The improved stochastic diffusive step Axjiff (<f,) equals 



Axdiff(^,) 



^4iff + ^-^diff + 2Axdiff 



(25) 



Here Axjiff is the rms diffusive step l |13l l at the old position and 
Axj.jj corresponds to the rms diffusive step evaluated at the two 
supporting positions x+ that were defined in relation M2V . 



Ax±ff = V2Z)(x±)Af . 



(26) 



The second term of the corrected diflFusive step J25b corrects for the 
'lopsidedness' of the random walk that results from the gradient in 
the diffusivity. 



4 A PREDICTOR-CORRECTOR SCHEME FOR SPATIAL 
TRANSPORT 

We have found that a scheme proposed by Kloeden & Platen ( 1 992) 
gives excellent results in the case of strong drift due to gradients in 
the diffusivity, where the CES fails to be sufficiently accurate. It is a 
second-order predictor corrector scheme, called the KPPC scheme 
in what follows. In particular this scheme improves the accuracy of 
the spatial transport of the particles, which leads to a better sam- 
pling of the acceleration rate w(x) in the shock compression. For 
the problem at hand this scheme takes the following form: 



Step 1: first supporting position value 

As a first step one calculates a first supporting position value x us- 
ing the Cauchy-Euler scheme: 



■+U{x, t) At + ^j2D(x)Kt^, 



(21) 



For simplicity we adopt a constant time step A/. The stochastic vari- 
able ^, 6 N(0, 1) is drawn from a normal distribution with zero 
mean and unit dispersion. Two additional supporting position val- 
ues are calculated. 



X + U{x , Af ± V2£>(x)Ar , 



(22) 



that correspond with the position of two hypothetical particles that 
experience ± the rms diffusive step. The stochastic variable ^, is 
now fixed at the value used in ( I21b . It does not change in any of the 
subsequent steps of the algorithm that are outlined below. 

Step 2: predicted position value 

As a next step one calculates the predicted position x at time f -1- Af 
as 



X = X + (/ Af + Axdiff (f,) . 



(23) 



The mean velocity U used in the advective + drift term is an aver- 
age velocity, defined by using the first supporting value x: 



= { [U{x) + U(x) ] 



Step 3: corrector step and final position 

One finally obtains the corrected (and final) position x^ = x(t + At) 
at time t + At in the following way: 



x,=x+^ [t/(x) + U(x) ] At + Axj,ff(^,) , 



(27) 



This last step uses the same diffusive step as in the predictor cy- 
cle but corrects the advective step using the predicted position 
X (see Eqn. I23t . This is a reasonable approach as (on average) 
Axdiff » Axadv. It was already noted by Marcowith & Kirk (1999) 
that a careful treatment of the advective step (including drift) is 
more important for the accuracy of the scheme than the treatment 
of the diffusive step. This scheme and the tests presented below 
bear that out. 

A few remarks about the implementation of this scheme are in 
order. 

First of all, this scheme is computationally about 6 times more 
expensive than the Cauchy-Euler scheme. An order-of-magnitude 
increase of the computational effort is typical when switching from 
an explicit, first-order scheme to a second-order accurate predictor- 
corrector scheme or the closely related Runge-Kutta type schemes. 

Secondly: for the term that corrects for diffusivity gradients 
to be effective one should not employ an often-used numerical 
approximation for the Wiener process that replaces the normal 
distribution for ^, by a two-point distribution of values, choosing 
^, = ±1, where the two possible signs are drawn randomly with 
equal probability = P_ = |. In that numerical approxima- 
tion for the Wiener process the second term in the expression ( 125 1 
vanishes identically, and much of the scheme's improved accuracy 
with respect to the Cauchy-Euler scheme is lost. In that respect one 
might expect that a symmetric three- value scheme for for exam- 
ple (in the notation [value | probability]) 



■V3| i 

D 



2 

°l3 



■V3|^ 
o 



(28) 



(24) 



works better. 
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5 TESTS OF THE ALGORITHM 
5.1 Basic assumptions 

We have tested the KPPC scheme as described here, comparing 
its performance to the performance of the simpler Cauchy-Euler 
scheme. For this test we use scaled (dimensionless) variables where 
the fluid velocity V is measured in units of the shock speed and 
position X along the shock normal is in units of the shock thickness 
Ls For clarity we keep in the equations even though Lj = 1 in 
the numerical implementation. The velocity V{x) is in the direction 
of positive x, given by 



V(x)-. 



r+l 



r-\ X 

2r- — ''^"^ z: 



(29) 



The velocity decreases with increasing x, from V(-oo) = Vi = 1 
to V(+oo) = V2 = 1 jr. This means that we work in the rest frame 
of the shock and measure the flow velocity in units of the shock 
velocity with respect of the upstream medium. Here r > 1 is the 
compression ratio of the shock transition in the sense that (for 
this one-dimensional steady flow) the conservation of mass implies 
pV = constant, with p the mass density. The density contrast be- 
tween the far upstream and far downstream state follows as 

p(oo) ^ P2 _ Y± _ 
p(-oo) pi V2 



(30) 



This velocity profile models the shock as a stationary and smooth 
transition, with a width (velocity gradient scale) Lj. To model a 
varying dilfusion coefficient we adopt a diffusion coefficient that 
varies with position x as 



D(x) = Di 



0-+ 1 



2a- 



^-1 

tanh — 

2a- \U 



(31) 



Here Di is a constant dimensionless diffusivity that is related to 
the physical diffusivity /)phys far ahead of the shock by D[ = 
Dp^yJLsV, with V, the shock velocity. The diffusivity decreases 
if one moves from upstream (x < 0) to downstream (x > 0) across 
the shock, with a ratio of asymptotic values equal to 



D(-oo) £>i 

= — = O" > 1 , 

£)(oo) D2 



(32) 



the kind of behavior one expects in astrophysical applications. The 
scale length for the variation of the diffusion coefficient in these 
units is of order Ld. For future use we define the quantity 



(33) 



Skadron, 1977, Bell, 1978; Blandford & Ostriker, 1978). In present 
notation, using the momentum p rather than y = ln(p/rac): 

dN _„ , „, 3 



N{X, p) : 



Ax d In p 



P 



q(s = 0) 



r-l 



(35) 



In our tests of the algorithm we have assumed a diffusion coeffi- 
cient that is independent of particle momentum, so this power-law 
behavior is valid uniformly across the grid, with only the concen- 
tration of test particles varying with position x. 



5.2 EfTect of finite shock thickness 

Since we assume a finite shock thickness, a situation typical of 
shocks obtained through numerical simulation, we need an ex- 
pression for the slope q for finite e. We use a perturbation analy- 
sis adapted from Drury (1983) and the closely related method of 
Schneider & Kirk (1987). The analysis presented below is valid 
when there is a small parameter e, in this case the ratio of the shock 
thickness and the cosmic ray diffusion length: 



« 1 



(36) 



Consider the steady-state transport equation (|3) reformulated in 
terms of the Vlasov distribution f{x , p) = dN/(d^x d^p). In our 
application we have f{x , p) = N(x , p)/4np-^. Assuming a one- 
dimensional steady flow in the x-direction the equation for f(x , p) 
with df/dt = reads: 



dx dx \ dx 



1 dV 
3 d7 



df 
'dp 



(37) 



Schneider & Kirk use a variant of the Ricatti transformation (e.g. 
Polyanis & Manzhirov, 2007, Ch. 12.2), which we slightly modify 
here. We introduce a dimensionless position variable X, 



X(x): 



r V(x')dx' 
Jo D(x') 



(38) 



and define (c.f. Schneider & Kirk, 1987) 



D 



G(x, p) = - 



df 
dx 



V 



dX 



f(x , p) f(x , p) 



(39) 



Schneider & Kirk also assume a power-law momentum depen- 
dence. 



f(x , p)oc p 



(40) 



This is essentially the Peclet number of the shock based on the 
cosmic ray diffusivity. In terms of this quantity one has 



Axd, 



Ax„, 



2a- 



1 L, 



(34) 



A sharp shock in the present context corresponds to e <c 1. 

The main test of the algorithm lies in its ability to reproduce 
the spectrum predicted by the analytical theory of DS A at a steady 
shock with e ^ I, cr ^ r and - Lj. In the limit of a infinitely 
thin shock with £ = (in physical terms: a shock thickness that is 
much smaller than the scattering mean free path of the accelerating 
cosmic rays) and in the absence of radiation losses the predicted 
shape of the spectrum is a power law in momentum, with an index 
q that depends only on the compression ratio r (e.g. Axford, Leer & 



which can only be strictly justified if the diffusion coefficient is 
independent of momentum, the case we consider here, and for 
momenta well above the injection momentum. Note that we have 
^ = ^ -I- 3. In that case G(x , p) is a function G(X) of position alone. 
Equation l l37b can be written as a non-linear ordinary differential 
equation for G(X): 

(41) 



dG q dV 
dx'^ 3dX 



G + 



This is the relation derived by Schneider & Kirk (1987), general- 
ized to the case of a position-dependent diffusion coefficient. The 
boundai-y conditions for G(X) alX = ±00 are 



G(-oo) = -V(-oo)= -Vi , G(+oo) = 0. 



(42) 
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The first condition assumes that there are no pre-existing particles 
far ahead of the shock so that asymptotically f(x , p) oc exp(X) 
for large negative X, where V(X) =^ Vi is approximately constant. 
The second condition, which states that the diffusive contribution 
oc dfldx to the flux vanishes asymptotically far behind the shock, 
ensures that the particle density remains finite as X — » +00. 

We now assume that condition ( 1361 ) is satisfied, meaning that 
the shock transition must be sufficiently sharp. In that case, the left- 
hand side of Eqn. ( I41t will be large in a region AX ~ e when com- 
pared with the two terms on the right-hand side. This behavior can 
be formalized by using e as a formal ordering parameter, replacing 
(Eg by 



1 AG q dV\ ^ G- 
- — + - — \ = G+ — 
e\AX 2, ax V 



We seek solutions of the form 

G{X) = Go(X) + £ Gi (X) + £- GiiX) + ■ 

and expand the slope as 
q = qo + s qi + s' qi + ■ ■ ■ 



(43) 



(44) 



(45) 



We can now solve ( 143 1 at each order of s, putting £ = 1 at the end 
of the calculation. At leading order (£"') one has 



dGo _^ 90 dV ^ Q 
AX 3 AX~ ' 

subject to the boundary conditions ( I42t for Go(X): 
Go(-oo) = -Vi and Go(+oo) = . 



(46) 



(47) 



Integrating i52\ from X = -00 lo X = +00 immediately yields a 
relation for qi : 



J (V2- V,): 



:£"dx(Go4 



Solving for qi, using ( I50t and ( I5U : 

Vi (Vt -V)(V- V2) 



qi = qo AX 



i: 



V (Vi - V2Y 



(54) 



(55) 



The function Gi (X) is 



G,(X) = ^(Vi-V) 



(56) 



AX 



, Vi V2(Vi-V')(V'-V2) 
V (V, - V2)- 



Here V = V{X'). Note that qi and Gi(X) vanish automatically 
if one uses the step function velocity profile of an infinitely thin 
shock, with V{X) = V, for X < and V(X) = for X>0. 



At order £ one finds the following equation for G2(Xy. 
AG2 h AV ^ / 2Go\ 



(57) 



subject to the boundary condition G2(-'») = G2(+o°) = 0. Inte- 
grating ( I57t from X = -00 lo X = +00 yields an equation for q2: 



The solution is elementary: 
Go(x) =j[V2- V(X) ] , 



(48) 



with V2 = V(+oo) the asymptotic velocity downstream. The zero- 
order slope qo can be found by integrating ( 146) from X = -00 to 
X = +00 and using the boundary conditions. Another elementary 
calculation gives 



^ (V2 -Vi) + Vi=0. 



(49) 



One finds (as expected) that qg is the slope associated with an in- 
finitely thin shock, where the velocity jumps from Vi to V2 < Vi at 
X = 0: 



qo = q(s = 0) 



3V, 



V, - 



(50) 



V, - V- 



-j dZGi(X)| 
-J dXGi(X)| 



2Go(X) 
V(X) 



(58) 



(Vi + V2) V(X)-2ViV2 \ 

v(x)(yi-V2) 



This procedure can be extended to higher order, but little is gained 
at the expense of increasingly complex mathematics. 



5.2.1 Specific examples 

We use two examples of immediate importance for a test of the 
numerical scheme advocated here. Table 1 gives the parameters as 
used in the numerical simulations presented below. 



This also ensures the the boundary condition atX = -00 is satisfied 
as i50\ substituted into ( 148 1 implies 



Go(X) = 



Vi (V(X) - V2) 



Vi - V2 
At next order (£") one has: 

^ + ii^ = G +^ 

AX 3 AX ° V ' 
subject to the boundary condition 
Gi(-oo) = Gi(+co) = 0. 



(51) 



(52) 



(53) 



Model 1: hyperbolic tangent velocity profile and constant 
diffusivity 

The first example is the case of a uniform diffusivity D{x) = D\, 
which formally corresponds to tr = 1 and = 00. This is the case 
where the CES is known to yield good results. This case has been 
treated before by Axford, Drury & Summers (1982), who show 
that for the hyperbolic tangent velocity profile ( 129) adopted here 
the cosmic ray transport equation can be solved analytically. They 
find an asymptotic slope of the momentum distribution equal to 



?o 1 + 



V2L, 
2D, 



(59) 
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Here qo = 3r/(r - 1). The perturbation expansion used here (and 
in a slightly different form by Drury (1983)) reproduces this result. 
The hyperbolic tangent velocity law M9\ implies 



dV 

dx 



(60) 



Substituting this into the generally valid expression ( I55K together 
with AX = V dx/Z)i, one finds: 



9i 



20i(Vi - Vi) 



This agrees with result ( 159) of Drury et al (1982). Using this in 
relation ( 156) one finds that Gi(X) = 0, which implies G„ = q„ = 
for n>2. Here the perturbation expansion breaks off at order £ and 
yields the exact asymptotic result, as noted before by Drury (1983). 
In terms of the compression ratio r = V1/V2, the diffusion length 
far upstream L^^ = Di/Vi and q = q - 3 one has: 



-^fl.^] = ^(l.f) 
r-l \ IDij r-\ \ 21 



(62) 



Model 2: hyperbolic tangent profile and constant diffusion length 

As a second example we consider the case of a constant diffusion 
length: 

_ D{x) _ Oi 
~ V(x) ~ V, ■ 



(63) 



This example is important as a test case of the predictor-corrector 
algorithm used here as it has a strong gradient in the diffusivity. 
Formally it corresponds to cr = r and = L, so that in the shock 



Axd, 



1 



2rB 



(64) 



which becomes large if £ « 1 for thin shocks. If one adopts the 
hyperbolic tangent profile (|29)/ll60t and uses the fact that 



dX: 



Ax 



(65) 



relation ( 155) yields: 



2idiff 



VI-V2 



AV 
VAx 



(66) 



2Ldiff \V1-V2 \V2 



Table 1. Model parameters 



Model r cr Li/L, Vi At/L, 

1 4 1 00 0.05 

2 4 4 1 0.05 



strong shock in a mono-atomic gas, one has r In r/2{r- 1) ^ 0.924. 
We have obtained the second-order term for the important case 
r = 4 through numerical integration of the integral in the expression 
for qo. We find: 



q(r = 4) 



1 + 0.924 — + 0.095 



1 + 0.924 £ -t- 0.095 £^ 



(68) 



Note that the end result in both cases is a series in £ = LslLm- 
It should be pointed out that the term oc £- is small in the second 
case, and vanishes completely in the first case. The smallness of the 
second-order coiTection to q seems to be a rather general property 
of this expansion in £ for reasonable velocity profiles. As a further 
example: in the case of a linear velocity profile, where 



V{x)-. 



V, 



Vi + V2 Vi- V2 



forx < -LJ2, 



X for -LJ2 <x< LJ2, 



for X > Ls/2, 



(69) 



and for a constant diffusion coefficient D = Di (a = I, = 00), 
the same procedure yields: 



3 /. £ (r+l)£- 

fl = fl-3^ 1 + - + — 

^ ^ r-l\ 6 360r 



For r = 4 this is 

£ 

q(r = 4) = I + - + 

^ 6 288 



(70) 



(71) 



These three examples suggest that the results obtained here for the 
slope q(E) are applicable even for £ 1 . The simulations presented 
in the next Section bear this out. 



The function Gi can be calculated, but the integral over Gi that 
determines the next order correction ^2 to the slope can not be 
expressed in elementary functions. Limiting ourselves to the first- 
order correction one has in terms of q = q - 3 and r = V1/V2: 



r- 1 



1 + 



r (In r) L, 
2(r- l)Ldif 



r- 1 



1 + 



r (In r) s 



(67) 



Here the steepening due to a finite shock width ~ Lj is more pro- 
nounced than in the previous cases as the diffusion length near the 
shock transition is smaller compared to the first case. For instance: 
in a shock with compression ratio r = 4, the value expected for a 



5.3 Numerical results 

The two figures below show the results of numerical simulations 
for the two cases listed in Table 1 . The approach is similar to the 
one employed by Achterberg and Kriills (1992): particles are in- 
jected close to the shock, and followed over a grid that extends 
from X = -10 to = -l-lO. Particles are detected as they cross 
the downstream boundary X^ra - V2Xmax/£^2 = 10, which acts 
as an absorber. The influence of a downstream absorbing bound- 
ary on the slope q decays as exp(-Xn,ax) with respect to unity, and 
is negligibly small for these parameters. Particle splitting is used 



8 Achterberg (ir Schure 



slope momentum distribution for D = constant 



slope momentum distribution for L.j„ = constant 



q = 1 + 0.5 e 

* = CES, * = KPPC 




0.01 



0.1 



Figure 1. Results for Model 1: the case of a hyperbolic tangent velocity 
profile and constant diffusivity. The solid curve gives the result j62t . the 
open stars are the results of the KPPC scheme and the solid stars the results 
obtained using the CES. The shock compression ratio equals r = 4, the case 
of a strong hydrodynamical shock. The parameter e = Ls/Ldiff varies from 
0.01 to unity. Note that the figure employs a logarithmic scale for e. 
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Figure 2. Results for Model 2: the case of a hyperbolic tangent velocity 
profile and constant diffusion length Ljiff = D(x)IV(x). The solid curve 
gives the result (68), the open stars are the results of the KPPC scheme and 
the solid stars the results obtained using the CES. The shock compression 
ratio equals r = 4, the case of a strong hydrodynamical shock. Again a log- 
linear scale is employed in this figure. Note that the vertical scale differs 
from the one used in Figure 1 . 



at intervals equidistant in log p to minimize the effects of Poisson 
noise at large momenta where fewer particles reside in the distri- 
bution. The spectra obtained in this way are strict power laws that 
extend over five decades (Ay ^ 11.5) in particle momentum. We 
use a fixed time step that corresponds to Vi A? = 0.05 L^, so that the 
advective step resolves the shock transition. In practice, good re- 
sults are obtained if Ax^dv ~ 0.1 Ls- The difi^usion coefficient varies 
from Di = \ to Di - 100, which corresponds to a diffusive step in 
the range Axjiff =^ 0.3-3. Note that in our implementation we have 
scaled the spatial coordinate x with the shock width so that = 1 . 

Figure 1 shows the results obtained for the case of a constant 
diffusion coefficient (Model 1). In this model there is no drift term. 
In this case the CES and the KPPC scheme give comparable results 
that closely follow the theoretical prediction up to £ = L^jL^^g = 
0.04. For smaller values of s (i.e. larger values of Di) the both 
schemes become inaccurate as they under-sample the acceleration 
rate ii)(x) in the shock transition. 

Figure 2 shows the results for Model 2, where L^jfj = 
D(x)/V(x) is kept constant so that Ls = Lj and cr = r = 4. This 
is the model with a large drift in the shock due to the gradient 
in the diffusivity. The results for the slope of the momentum dis- 
tribution obtained using the CES (the filled stars) deviate signifi- 



cantly from the slope obtained using the perturbation expansion if 
£ = LJL^iff = Lj/Ljiff < 0.1. In contrast, the KPPC scheme gives 
significantly better results that are usable up to £ ^ 0.02. The error 
in the value of the slope q retumed by the KPPC scheme is typically 
3 times smaller than the error produced by the CES. The deviation 
from the analytical result becomes significant if the magnitude of 
the drift term in the statistically sharp spatial step Ax, becomes of 
the same order as the shock thickness: 



do 

dx 



A/>L, 



(72) 



Making the estimate |dZ)/dx| ^ |A£)|/Ls as it applies to typical sit- 
uations where the diffusion coefficient jumps by an amount AD 
across the shock, the accuracy of the KPPC scheme is lost if 



Ax<i„f, |AD| ViAt 



> 1 . 



(73) 



For instance: in the results shown in Figure 2 the deviation in the 
slope returned by the KPPC scheme becomes large when lAx^rif,! =^ 

In our test of the KPPD scheme we have assumed that the 
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diffusion coefficient is independent of momentum. In many astro- 
physical applications one expects the mean free path to increase 
with momentum so that the diffusion coefficient scales as Z) oc p". 
For instance: one often assumes Bohin Diffusion with a mean-free- 
path equal to the gyro radius, A^f^ = ^ pc/qB where q is the 
particle charge. For relativistic particles (v ~ c) this implies D oc p. 
The KPPD scheme should be able to give reliable result in this case 
also, as long as the time steps are such that lAxdriftl ~ Lj. The scal- 
ing D oc p" implies Axj^ft oc p", so it may be necessary, depending 
on the dynamic range in p, to employ the scheme with smaller time 
steps for the high-energy particles in the population than for the 
low-energy particles. 



6 DISCUSSION AND CONCLUSIONS 

We have argued that the simulation of diffusive shock accelera- 
tion of cosmic rays through the numerical integration of stochastic 
differential equations needs a more accurate, second-order scheme 
when large gradients in the cosmic ray diffusivity arise. This situ- 
ation is astrophysically relevant for oblique shocks or for shocks 
where the cosmic rays generate strong magnetic turbulence in 
the shock vicinity. The Kloeden-Platen predictor-corrector scheme 
proposed here is such a scheme. We have demonstrated using sim- 
ple simulations that the Kloeden-Platen predictor-corrector scheme 
is significantly more accurate for small shock widths coupled with a 
strong gradient in the cosmic ray diffusivity. For large shock widths 
(Peclet numbers larger than ~ 0.25), or when the gradient in the 
cosmic ray diffusivity is small, the two schemes produce compara- 
ble results in terms of the accuracy of the momentum distribution 
that is obtained for the simulated particles. 

Given the fact that the Kloeden-Platen predictor-corrector 
scheme is computationally about six times more expensive than the 
Cauchy-Euler scheme, one might consider implementing a hybrid 
approach where one switches between the simpler Cauchy-Euler 
scheme and the Kloeden-Platen predictor-corrector scheme with 
the switch based on the value of lAxdrift/Lsl, the magnitude ratio 
of the drift term and the shock thickness in the stochastic differen- 
tial equation 1 11 It and the shock width. The results obtained here 
suggest that the KPPC scheme is needed for sufficient accuracy 
whenever |AX;j,.ift| > 4Ax.^^, close to the shock. 

We have checked whether results with similar accuracy can 
be achieved with the Cauchy-Euler scheme, simply by reducing 
the time step until the computational expense is similar to that of 
the KPPC scheme with the larger time step. This turns out not 
to be the case. As an example we consider Model 2 for the case 
e = 0.04. The analytical estimate for the slope is qth = 1.037. 
We ran both schemes with ViAt = 0.1, 0.05, 0.025 and 0.0125 
in units where Ls = 1, thus halving the time step each time. The 
KPPC scheme consistently returns (within errors due to Poisson 
noise) a slope ^kppc = 1.035, quite close to the (approximate) the- 
oretical result. Table 2 gives the corresponding result for the slope 
obtained with the CES, ^ces- The slope still has a sizable error 
even for the smallest time step, 8 times smaller than the largest step 
where the KPPC scheme already performs satisfactorily. The slope 
Ices, although decreasing as the time step gets smaller, remains 
consistently too large. We conclude that, for given computational 
expense, the KPPC scheme is still superior. 

We note in passing that the Kloeden-Platen predictor-corrector 
scheme may also be useful in other numerical applications of 
stochastic differential equations, such as the solution of the equa- 
tions that describe stochastic particle acceleration by waves (Fermi- 



Table 2. Performance Cauchy-Euler scheme 



Advective step V, Af/Ls 0.1 0.05 0.025 0.0125 
Slope ?CES 1-222 1.150 1.098 1.077 



The slope q of the simulated momentum distribution retumed by the 
Cauchy-Euler scheme for Model 2 with e = 0.04 for different time steps. 
The theoretical slope equals q^^ = 1.037. 

II acceleration). Here the equation for the evolution of the momen- 
tum distribution of the accelerating particles contains a momentum 
diffusion term that can be written in the form (e.g. Melrose, 1980) 



df\ _ _}_d_l , d£ 



(74) 



and a large drift term is unavoidable. The relevant drift velocity (in 
this case corresponding to the mean momentum gain) is 

(I) =^1{P'^^)- (75) 
\ P- dp ^ " 

In these expressions Dp is the momentum diffusion coefficient. 
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